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Abstract. The main goal of this paper is to present results 
of comparison study for the level set and direct Lagrangian 
methods for computing evolution of the Willmore flow of 
embedded planar curves. To perform such a study we con- 
struct new numerical approximation schemes for both La- 
grangian as well as level set methods based on semi-implicit 
in time and finite/complementary volume in space discretiza- 
tions. The Lagrangian scheme is stabilized in tangential di- 
rection by the asymptotically uniform grid point redistribu- 
tion. Both methods are experimentally second order accurate. 
Moreover, we show precise coincidence of both approaches 
in case of various elastic curve evolutions provided that solv- 
ing the linear systems in semi-implicit level set method is 
done in a precise way, redistancing is performed occasion- 
ally and the influence of boundary conditions on the level set 
function is eliminated. 



which is a critical point (minimizer) for the elastic energy 
functional 



(1) 



The first comprehensive study of analytical properties of non- 
closed planar curves that are minimizers to ([T]i goes back 
to Leonhard Euler who presented their characterization and 
classification in the pioneering work Additamentum I (De 
Curvis Elasticae) contained in his Opera Omnia |11|. Since 
then much effort has been spent to analyze and provide com- 
plete characterization of both minimizers to ([T]i as well as so- 
lutions corresponding to the gradient flow associated with the 
elastic energy functional ([T]i. It is well known from Euler's 
work that the flow of planar curves with the normal velocity 
given by 



/3 



1 



(2) 
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1 Introduction 

In the past years, elastic curves, the Willmore functional and 
the corresponding gradient flow (the Willmore flow) attracted 
a lot of attention from both theoretical as well as computa- 
tional point of view. Following Daniel Bernoulli's model of 
an elastic rod, a classical elastica is a curve F in the plane 
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is a gradient flow for the elastic energy functional E{r) (see 
e.g. 17,10]). Such fourth order flows of closed curves and its 
3D analogies appear in various physical and computer vision 
applications dealing with a motion of phase interfaces or with 
an image and surface reconstructions l5] fT0ll22lfT5]l 6ll4l l25l . 

We remind ourselves that the so-called surface diffusion 
problems (see e.g. |2,20|) are described by nonstationary 4*'' 
order intrinsic partial differential equations. Similarly, a nu- 
merical solution to the Willmore flow, either in direct (La- 
grangian) or level set (Eulerian) formulation, is a nontrivial 
problem and leads to a solution of fourth order in space non- 
linear evolution PDEs Convergence of a semidiscrete time 
continuous finite element discretization in the case when the 
evolved surface is a graph has been proved by Dziuk and 
Deckelnick in |8|. First numerical study based on the finite 
element method for the Willmore flow in Lagrangian formu- 
lation was presented in lITOl and for the level set formulation 
in |9|. Finite difference discretization has been analyzed in 
H.,21J. Tangential stabilization of Lagrangian approach for 
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solving fourth order elastic curve flows in case of surface 
diffusion was first introduced in |20|. Then a parametric fi- 
nite element method was tangentially stabilized in |3|. Al- 
though the Lagrangian methods are fast and robust (when in- 
corporating a suitable tangential velocity) they cannot handle 
topological changes for which the level set methods are pre- 
ferred II22I I9I . However, a careful and systematic comparison 
of nontrivial examples of direct and level set approaches for 
fourth order curve evolution problems is still missing. The 
goal of this paper is to provide such a comparison study, and, 
moreover to derive new numerical schemes based on the fi- 
nite/complementary volume strategies for both Lagrangian 
and level set formulations of the Willmore flow. 

First, we present a tangentially stabilized Lagrangian met- 
hod based on a solution to the curvature, local length and 
position vector equations accompanied by the asymptotically 
uniform tangential redistribution of numerical grid points. We 
show experimentally that the method is second order accu- 
rate. We apply this method to various examples of evolution 
of planar embedded curves. Stabilization by the tangential ve- 
locity allows us to use reasonable large computational time 
steps and prevent formation of various instabilities like merg- 
ing of evolving curve representing grid points or swallow 
tails, which are typical disadvantages of the direct methods. 

Then we develop new semi-implicit complementary vol- 
ume scheme for solving level set formulation of the Willmore 
flow. It is again second order accurate. Due to a finite vol- 
ume character of discretization it has a potential to be nat- 
urally connected with finite volume schemes for advective 
level set equations II121I13I and thus to be used in various 
models where the fourth order terms serve as a curve motion 
regularization arising, e.g., in image segmentations [25 1. 

The outline of the paper is as follows. In section 2. 1 we re- 
call a closed governing system of equations for the curvature, 
local length and position vector describing evolution of plane 
curves satisfying (|2]l in Lagrangian formulation and describe 
the main idea of asymptotically uniform tangential redistribu- 
tion. Section 2.2 focuses on the brief derivation of the govern- 
ing equation representing the evolution of zero level set sat- 
isfying the geometric equation (|2]l. In section 3.1 we present 
our Lagrangian numerical approximation scheme and, in sec- 
tion 3.2, approximation of the level set equation for the Will- 
more flow. Section 4 is devoted to study of the experimental 
order of convergence for both methods and to comparison of 
both methods in various elastic curve evolution examples. 

2 Governing equations 

2.1 Direct Lagrangian method 

Henceforth we shall parameterize an embedded regular plane 
curve r hy a. smooth function x : ^ R^, i.e. F = 
Image(a;) :— G S^} such that the local length el- 

ement g = \dux\ > is everywhere positive. Taking into ac- 
count the periodic boundary conditions at ?i = 0, 1 we shall 
hereafter identify with the interval [0,1]. The unit arc- 
length parameterization will be denoted by s, s,o As — g du. 
The tangent vector T and the signed curvature k of F satisfy 
T — dgX = dux/g, k ~ dgX A d'^x — duX A d^x/g^. More- 
over, we choose the unit inward normal vector N such that 



T A N = 1 where a A b is the determinant of the 2x2 matrix 
with column vectors a, b. By f we denote the tangent angle 
to F, i.e. T = (cos i^, sinz^)-^. Now it follows from Frenet's 
formulae that ^^T = fcN, S^N = -kT and dsiy = k. Notice 
that the curvature k is positive for convex closed curves in our 
convention of picking of normal and tangent vector orienta- 
tion. 

Let a regular smooth initial curve Fq = Image (xq) be 
given. According to ifTTll . an evolving family of planar curves 
Ft = Image(a::(., i)), t € [0,T), satisfying Q can be repre- 
sented by a solution to the following system of PDEs 

dtk = d^f3 + adsk + k^P, (3) 

dtg = -gkP + gdsa , (4) 

dtx = /3N + aT , (5) 

subject to initial conditions k{., 0) — ko , g{., 0) = go , and 
x{.,0) = xo{.),. We impose periodic boundary conditions 
at u = 0,1. Having recalled the general form of governing 
equations we are able to calculate the time derivative of the 
elastic energy functional 

2—E{Ft)^—[ k^gdu= [ 2kdtk ^ k^ (3 + k^dsods. 
dt Jq J 

Since /^^ kdgf3 = Jp^ Pd^k and Jp^ ds{ak^) we obtain 
the following equation 

^^E{Ft)^ I^{d',k+^k^)f3ds. (6) 

It enables us to conclude that the evolution of Ft with the 
normal velocity f3 — —d^k — ^k^ is a gradient flow (the 
Willmore flow) for the Willmore elastic energy functional E. 

Notice that the tangential velocity a is a free parameter in 
(O-Q and it may depend on other quantities like e.g. the cur- 
vature, normal velocity and/or local length element in various 
ways including local or nonlocal dependences, cf. [14..16.iT7l 
[T8l[T9l . In this paper we make use of the so-called asymp- 
totically uniform tangential redistribution derived in 111811191 
which is the most natural for the Willmore flow since an ini- 
tial shape is approaching evolution of expanding circles. Let 
us denote L = Lt the total length of a curve Ft. It follows 
from analysis of the tangential velocity made in 111811191 that 

g{u,t) -, 

— > 1 as t ^ Tmax uniformly w.r. tou G S 

Lt 

provided that the tangential velocity is a solution to a non- 
local equation 

d^a = kf3- {k(3) r + [L/g - 1) a(0, t) = Q. (7) 

Here w > is a given positive constant and {.)r is an aver- 
aging operator over a curve F, i.e. (fc/J) r — Jp kf3ds. It is 
clear that redistribution of grid points along a curve becomes 
uniform as t approaches the maximal time of existence Tmax- 
In the case of a Willmore flow the time horizon is infinite (i.e. 
Tmax = +00) for planar Jordan curves and Tnax can be fi- 
nite for some selfintersecting immersed curves in the plane. 
Furthermore, inserting a computed from (|7]l into ©-(ISll and 
making use of the identity adgk = ds{ak) — kdga then the 
curvature and local length equations can be rewritten as fol- 
lows 

dtk = dlP + ds{ak) + k{k/3) r + (1 - L/g) kuj , (8) 
dtg = -g{k(3)r + {L~g)Lj. (9) 



COMPARISON STUDY FOR WILLMORE FLOW OF PLANAR CURVES 



3 



In other words, the strong "point-wise" influence of the term 
kf3 in ^ and (|4]i has been softened by the "averaged" term 
{k/3) r in (O and (|9]l. As a consequence, this important prop- 
erty of asymptotically uniform tangential velocity enables us 
to construct an efficient and stable numerical scheme prevent- 
ing fast local decrease of local lengths (merging of numerical 
grid points) as well as forming various further numerical in- 
stabilities related to high local curvature. Since 

d^x = d^T = dl{kN)dlkN + 25,fc9«N + kd^N 
= affcN - 2{dsk)kT - kds(kT) 
= affcN - 3fc(9,/c)T - k^dsT 

= affcN - '^ds{k^)dsx - k^dlx 

and d^x = kN we have 

i-d^k ~ h')N = -dtx - h^dlx - ^dsik^)d,x 

= -dtx-^ds{k^d,x). 

Thus the governing system of equations (|3ll-(|5]l for the Will- 
more flow dU with tangential redistribution can be written as 
follows: 

dtk ^ -d^k - Idlik^) + ds{ak) + k{kp - d^a), (10) 



dti] = — fc/3 + dga, T] = ln(5), 
3 

dtx — —dgX — —ds{k'^dsx) + adgX 



(11) 
(12) 



where the tangential velocity a is the unique solution to equa- 
tion 



2.2 Level set method 

In the level set method the evolving family of planar curves 
^t, i > 0, is represented by the zero level set of the so-called 
shape function w : /2 x [0, T] ^ R where C is a simply 
connected domain containing the whole family of evolving 
curves ri , i 6 [0, T]. Assuming zero is the regular value of 
the mapping t), i.e. I Vu(a;, i)| ^ {)for:u{x,t) = Owe can 
express the unit inward normal vector and signed curvature 
as: N — Vu/|V7i| and k = — div (Vu/|Vu|). Let us denote 
the following auxiliary functions: 



= div 



Q = |Vu|, w = QH. 



Then dgk = —VH.dgX — — ViJ.T and, by Frenet's for- 
mula, d^k = -kS/H.N - T'V^i/T. Differentiating the 
equation u{x{s,t),t) = with respect to time we obtain 
dtu + Wu.dtx = 0. Since the normal velocity of x is (3 = 
dtX-N we obtain j^^dtu = —(3. Inserting expressions for 

d^k and fc = — _ff we obtain 

—dtU -div -TT-Vu \ -H'^ - TW^HT . 

Q \ Q J 

Here T — (— rt2,rti) where N = (rii,7i2), i.e. T is the vec- 
tor N rotated by —tt/2. Straightforward calculations show 
that the right hand side of the above equation can be rewritten 
in the divergent form. The resulting system of two equations 



governing the evolution of the shape function has been de- 
rived by Droske and Rumpf in ||9l and it reads as follows:: 



1 

dtu = -Qdiv ( EVw - T^753 Vm 



Qdiv 



(13) 
(14) 



I - ^ ^ I IS a projec- 



where the 2 x 2 matrix ]E — ^ (l\ 

tion into a tangential space of the curve representing the zero 
level set of u. System of equations ( fT3] - [l4l) is subject to the 
initial condition 

u{x, 0) = u°{x) , X e n 

and clamped boundary conditions at dH, i.e. u{x, t) — 0, 
d,yu{x,t) = 0, X € df2. The initial function u° is a signed 
distance function, i.e. — dist(a;, F^), x E il. 



3 Numerical approximation scliemes 

3.1 Numerical approximation of the Lagrangian method 

Our numerical approximation of an evolved curve is repre- 
sented by discrete plane points xj where the index i — 1, ...,n, 
denotes space discretization and the index j = 0, m, stands 
for a discrete time stepping. Due to periodic boundary con- 
ditions we use additional values x■'_^ — xi_T, xi — xl. 



a;jj_|_j^ — x\, a;^_|_2 = xl,. If we take a uniform division of 



the time interval [0, T] with a time step r = ^ and a uni- 
form division of the fixed parameterization interval [0, 1] with 
a step h — 1/n, a point x^ corresponds to x{ih,jT). The 
systems of difference equations corresponding to (|7]i, ( fTol i - 
( fT2] i will be solved for discrete quantities aj, rjf, , kj, xj, 
i = l,...,n, j = 1, m, representing approximations of 
the unknowns a, rj, gh, k and x, respectively. Here al rep- 
resents the tangential velocity of a flowing node xj , and rjf , 
rl « \xj—xj_i\ and kf represent piecewise constant approx- 
imations of the corresponding quantities in the so-called flow- 



ing finite volume 



In order to derive new position 



we use corresponding flowing dual volumes 



-1' ■ 



where x-i 



with approximate lengths q^- 



xl_i\. Our computational method is simple and natural. At 
the j-th discrete time step, we first find values of the tangen- 
tial velocity aj by discretization of (|7J. Then the values of r/j 
are computed and used for updating local lengths by dis- 
cretizing equations (fTTT i. Using computed local lengths, the 
intrinsic derivatives are approximated in JTOl i, and (fTZt . and 
pentadiagonal systems with periodic boundary conditions are 
constructed and solved for new discrete curvatures kj and po- 
sition vectors x^ . 

In order to discretize O we integrate it over flowing finite 
volume [xi-i,Xi] to obtain 



dgads 



kl3- {k(3)r+uj{L/g-l) ds . 
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Hence 

Ui - Ui-i = ri{ki(3i - {k(3)r) + {L/n~ n) 

where ao — 0. Taking discrete time stepping in the previous 
relation we obtain following expression for discrete values of 
the tangential velocity: 



Now replacing the time derivative by the backward differ- 
ence, derivative terms on boundaries of dual volume by fi- 
nite differences and Xi by the average of grid points in the 
last term, we obtain two tridiagonal systems for updating the 
position vector: 



(17) 



where, for i—\ 



, . . . , I c/, 




1=1 1=1 



and — 0, i.e. the point is moved in the normal direction. 

Now, a similar approximation methodology is applied for 
equation (fTTT i. Thus 



7 7^1 



for z = 1 , . . . , 71. It leads to the update formula for local lengths: 



exp(?7,^ 



I = 1, 



.,n, 



subject to periodic boundary conditions r^_-^ 

fin ffi+i = ^ij ~ ^2- local lengths are used for 

approximation of intrinsic derivatives in the curvature equa- 
tion ([lOll. We obtain /^"';_^ dtk ds = J^l^ -d^^k-\dl{k^) + 

ds{ak) ds + J^' ^ k{k[3 — dgo) ds. Hence 

+ ["^IxLi + hinkiPi - {ai - Qi-i)) (15) 



and taking semi-implicit time stepping, i.e. replacing time 
derivative by backward difference and treating linear terms 
at the current time level j while the nonlinear terms at the 
level j — I, and approximating derivative terms on the bound- 
aries of flowing finite volumes by finite differences we obtain 
following pentadiagonal system with periodic boundary con- 
ditions for new discrete values of the curvature: 



aiiki, + lr>ki_, + c>ki + d^kl, + 4kl^ = ff (16) 

for i = 1 , . . . , 71, subject to periodic boundary conditions k-'_^ = 

K^K^ K+i = = For completeness, a 

detailed description of the system coefficients is given in Ap- 
pendix. 

Finally we discretize equation (fT2] i by integrating in a 
dual volume Xi] to get 



dtxds 



dxi 



xi+i 2 

—d^x — -ds{k^dsx) + adgX ds , 



dt 



'dgX — -k^dgX 



-1 Xi+i 



+ a^{x.^+l - Xj) . 



1, 71, subject to periodic boundary conditions xi 



''n-l' -^0 



''m -^n+l 



li ■^n+2 



The exact form of 



coefficients A, B, C, V, £ can be found in Appendix. 

The initial quantities for the algorithm are computed from 
discrete representation of the initial curve xo- The reader is re- 
ferred to 1 19| for further details. Every pentadiagonal system 
is solved by Gauss-Seidel iterates. We stop the Gauss-Seidel 
iteration procedure if a difference of subsequent iterates in 
maximum norm is less than the prescribed tolerance lO^^*'. 



3.2 Numerical approximation of the level set method 

Concerning approximation of the level set equation ( fTSl l we 
consider rectangular domain = (01,02) x (&i,62) and we 
assume an equidistant spatial step h in both directions. We 
define a regular mesh ojh consisting of grid points Xij = 
[ai + ih,bi + jh] fori — 0,...,Ni,j — 0,...,N2, where 
fli + Nih = a2 and bi + ~ 62- Without loss of gen- 
erality we shall assume ai = 61 = 0. The corresponding 
dual mesh V is given as the union of the finite volumes Vij of 
the form {{i ^ ^) h, {i + ^) h) x {{j - ^) h, {j + ^) h) for 



1 = 0, Ni,j — 0, A'2. The projection of a solution at 
Xij is defined as Uij — u{xij ). Similarly as in the Lagrangian 
method we take a uniform division of the time interval [0, T] 
with a time step t = Let us consider an element Vij of 
the dual mesh V. Integrating (fT3Ti-(fT4ll over Vij and applying 
the Stokes theorem we obtain 



A ,- Qdt Jg 



1 du , . 

2 0^5^-^^^-'^^'^^' 



-^dx — 



-—da 

avu Q 9" 



(18) 



(19) 



where is the outer normal of the boundary dVij. 

We start with approximation of the term Q on Vij. For 
r,s € {—1, 1}, li'l + |s| = 1, we define the linear operator 
as follows: 



V'^Uij = i [r{ui+r,j 



Uij),U. 



r,l 



'1,S 



j+s 



where u^^ is the average of ity defined as: 



1 



^i+r.j 



For a fixed regularization parameter < e <C 1 we define 

1 

4 



Q 



0- 



E 

\r\ + \s\ = 
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Let E^J'" = (^^fef.'/J ) be the 2 x 2 projection matrix: Table 1. EOC for the Lagrangian scheme in ,p = 2, oo norms. 



-\rs\n 



-jrs;n 
■"ij 



Now we are able to derive a discretization of ( fTSb 



-jn— 1 



|r| + |s| = l 



^n — 1 



^;i(E^^"-V-<,zA..) (20) 



where 



|r| + |sj = l 



i-\-r.j-\-s J ' 



(21) 



r, s) for I 



and Vrs is the unit outer normal vector, — 
\s\ = 1. In order to approximate w" and w"^^ 
its boundary dVij we have used expression ( fT9] l to obtain 



on and on 



1 



O 

1 



(22) 



Since (l20b contains a new time level expressed through 
new time level of the solution (see (|22] |) as well as the pre- 
vious time level w"^^, the resulting discrete level set scheme 
is semi-implicit in time. After some calculations it can be 
written as twenty one points scheme of the form 



(r%s)eO 



(23) 



where O = {(r, s), -2 < r, s < 2, \r\ + \s\ < 4} and i ^ 
2, A^i — 2,j = 2, A'2 — 2. For the remaining the 
values of ufj are linearly extrapolated. The coefficients of the 
above system can be found in Appendix. 

System (|23] ) is solved by the iterative GMRES algorithm 
with ILUT (ILU with threshold) preconditioning or by the 
complete LU decomposition as a direct solver (c.f. ll23l ). The 
time step t is chosen to be proportional to and the regular- 
ization parameter e can chosen either as a function of /i or e 
can be prescribed as a small fixed constant. As an initial con- 
dition we choose a signed distance function d{x) to the initial 
curve. At prescribed redistancing time steps we perform re- 
distancing of the level set solution back to the signed distance 
using the fast sweeping method (see 1241 for details). As an 
alternative to the semi-implicit scheme (l20b we may also con- 
sider its explicit version, i.e. all the terms on the right hand 
side of ( l20b are considered at the time step n — 1, c.f. 1 1 ,21 1 
for other similar explicit schemes. In this case we avoid the 
singularities of the signed distance function in its local ex- 
trema and the initial condition has a "phase-field" like shape 

u°{x) = 6sgn{d{x)){l - exp{~\d{x)/S\), 

where (5 is a parameter describing the width of the region 
where u° changes from —5 to +5. 



Error \ r 

mc 


0.1 


0.05 


0.025 


0.0125 




0.04301 


0.01089 


0.00271 


0.00067 


EOC 




1.982 


2.005 


2.003 


p — 00 


0.03402 


0.00886 


0.00223 


0.00056 


EOC 




1.940 


1.986 


1.988 



Table 2. EOC for the level set scheme in L'',p = 2, 00 norms. 



Error \ ^ 
F.OC 


0.4 


0.2 


0.1 


0.05 


p = 2 


0.21497 


0.06585 


0.01699 


0.00400 


EOC 




1.707 


1.954 


2.086 


p — 00 


0.71190 


0.12286 


0.03780 


0.00973 


EOC 




2.534 


1.700 


1.957 



4 Discussion on numerical experiments 

4.1 Experimental order of convergence for the methods 

Let an initial curve be a circle with radius rp. Since for the 
circle we have k = ^ then it follows from ^ that r{t) — 

ir(t)~'^. Hence r(t) — (2t + ro) 3. Using this simple analyt- 
ical solution we can compute experimental order of conver- 
gence for both schemes. Without loss of generaUty we choose 
ro = 1. 

In the case of the Lagrangian scheme we approximate 
the initial unit circle subsequently by n = 10, 20, 40 and 80 
nodes with h = 1/n. The final time was set up to be T — 2.56 
and time step was chosen to be proportional to /i^, i.e. r = 
h^. Table [U shows errors and experimental order of conver- 
gence (EOC) of the scheme in LP = LP((0, T),Lp{S^)) = 
LP{S^ X (O, T)) for p — 2,00. In the level set approxima- 
tion we solve the problem in domain i7 — (—2,2) x (—2,2). 
The domain i7 was splitted subsequently into n x n finite 
volumes for n — 10, 20, 40, 80 with h = 1/n. The final time 
was chosen as T = 0.5 and again t = h^. The regularization 
parameter e was refined proportionally to the grid refinement 
using the rule = 2h. Finally, the redistancing period was 
Tredist = 0.25/i. Errors in LP,p — 2,oo norms are presented 
in Table 121 



4.2 Comparison of the Lagrangian and the level set 
evolutions 

In this section we compare the numerical results obtained by 
our Lagrangian and the level set approaches on various rep- 
resentative examples. In the case of Lagrangian scheme we 
approximate an evolving curve by 100 grid nodes in all ex- 
periments to follow. In the case of the level set method we 
hereafter split domain i7 into 100 x 100 finite volumes. 

In Fig. [U we present comparison of both methods for the 
case of evolution of an initial circle with the radius vq — 1. 
The time horizon T — 0.5. By cross marks we depict ap- 
proximation by the Lagrangian direct scheme where the evo- 
lution was computed using the time step r = 0.002 and 
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no tangential redistribution (a ~ 0). The evolution of the 
level set function was computed in the spatial domain H = 
(-2, 2) X (-2, 2} with the time step r = 0.002 and the 
smoothing parameter e — 0.001. We did not provide redis- 
tancing in this case in order to show deformation of an initial 
distance function to final shape of the level set function (see 
Fig.[T]bottom). 

In Fig.|2]we show evolution of an initial ellipse with half- 
axes 1 and 2. It asymptotically approaches a circle. We stop 
computations at the time horizon T — 2.56. In the case of 
the Lagrangian approach we pick r = 0.002 and the tan- 
gential redistribution parameter a; = 1. We computed evo- 
lution of the level set function in the computational domain 
i7 = (—4,4) X (—4,4). We chose the smoothing parameter 
e = 0.001 andwedidredistancingjustonceatTredist = T/2. 
Both, the initial and final level set functions are depicted in 

Figm 

In Fig. [3](top) the initial condition is a non-convex curve 
given by 




where < u < 1. The time evolution of such a non-convex 
initial curve was stopped at the time T — 0.01. In the La- 
grangian approach we picked r = 10^'', and the tangential 
redistribution parameter was u; = 1. The level set function 
was computed in the domain f2 = (—2, 2) x (—2, 2), with 
r = 2.5 • 10-5, e = 0.001 and Tredist = 0.001. Again, a 
comparison of the zero level set and the initial curve evolved 
by Lagrangian method show compatibility of both methods 
in the common time interval. In this example the Willmore 
flow quickly changes the shape of evolving curves from non- 
convex to a circular one. We show several time steps of the 
curve evolution in Fig. [3] 

In Fig. |4]we present evolution with the initial curve hav- 
ing sharp corners (see also a detailed close-up in Fig.|5]l. Al- 
though the initial curve (square) is convex, the evolved curve 
need not be convex for small times. Concerning numerical 
parameters, we chose n ~ 100 spatial nodes, r = 10"^ and 
the tangential redistribution parameter u — 1 (asymptotically 
uniform redistribution) in the direct Lagrangian method. As 
for the level-set method we took h = 0.04, r = 2.5 ■ 10^^, 
e = 10^5 and Tredist = 0.01. 

Fig.|6]shows comparison of the methods for another non- 
convex curve with very sharp corners. We chose the same nu- 
merical parameters as in the previous example for both the 
Lagrangian as well as level set methods. Also in this exam- 
ple one can observe satisfactory coincidence of numerically 
computed curves by both methods. 

Finally, in Fig. |7] we present an example illustrating a 
topological change. It has been computed by the level set 
method only because the direct method is unable to handle 
topological changes like pinching and splitting of curves. The 
initial zero level set consists of two almost touching curves - 
the inner curve being a circle and the outer curve being an 
ellipse with a shorter axis just slightly larger than the radius 
of the inner circle. We then let evolve this configuration by 
the level set equation. The phenomenon of pinching and sub- 
sequent splitting of the evolved curves can be observed in 
this example. Such a behavior can be observed in the mean 
curvature driven evolution of a dumb-bell initial surface in 




2 -2 

d) 

Fig. 1. a) A circle as an initial condition; c) the same initial con- 
dition for the level set approximation; b) a circle computed at the 
time T — 0.5 by both methods (cross marks correspond to the La- 
grangian method); d) the level set function at the time T — 0.5. 

3D where the Grayson theorem does not hold. To our best 
knowledge, there is no analytical proof of pinching-splitting 
phenomenon in the case of a Willmore flow of planar curves. 
Notice that this numerical result has been obtained only by 
using very small time steps, in the range 10^^^ — lO^^''. Since 
for such small time steps we do not increase efficiency by 
using the semi-implicit scheme we use here its explicit ver- 
sion with the Runge-Kutta-Merson fourth order adaptive time 
solver As further parameters we used h — 0.0166, e = lO^'^ 
and Tredist = 10^5. 
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a) 



b) 



a) 



b) 




c) 




c) d) 

Fig. 3. An initial condition a) given by i24l and curves computed 
by both methods at b): t = 0.001, c): t = 0.005 and d): t = 0.01. 




d) 

Fig. 2. a) An initial ellipse with half-axes ratio 1:2; c) the same 
initial condition for the level set approximation; b) the approximate 
circle computed at the time T — 2.56 by both methods (cross marks 
correspond to the Lagrangian method); d). the level set function at 
the time T = 2.56. 
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6 Appendix 

6.1 Coefficients of the Lagrangean systems 

The coefficients for the curvature system (fTSI l are as follows: 
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Using a similar strategy for approximation of the third order 
derivatives of position vector on boundaries of flowing dual 
volume we can write coefficients of (fTTI i: 
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6.2 Coefficients of the level set system 

The coefficients of the 21 -diagonal system matrix ( l23T l 
are given by: 



ill I ' 



{qlfrl^i rl{qlf rfqfqf_^ rj{qj_^)^ 

—-—— + — - ki Pi + — — , 

{ql-iYrl_^ ^ 2 2 



yl"" = 1 

11 



/l4 51 

r,se{-l,l} 



- rs;7i — 1 



2 q: 



,rs;n — 1 
ij 



+ 



E 



-\n~ 1 



\n—l 



l+|s|,2— |r|;2j I ^ — r, — s;n— 1 



q: 



4 l+|s|,l+|r|;ij I Q 



' 2+s j'+r 
-r, — s;n — 1 

i-\-s.j-\-r 



^ i — s.j—r 

Qsr;n — 1 
i — s.j — r 



where we used following approximation of third order deriva- 
tive terms on boundaries of flowing finite volume in (flSl i: 
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for I r I + I s I = 1 . Here we have denoted by Q* the harmonic 
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Next, for |r| 
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